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We provide an algorithm based on weighted-ensemble (WE) methods, to accurately sample systems 
at steady state. Applying our method to different one- and two-dimensional models, we succeed to 
calculate steady state probabilities of order 10 -300 and reproduce Arrhenius law for rates of order 
10 -280 . Special attention is payed to the simulation of non-potential systems where no detailed 
balance assumption exists. For this large class of stochastic systems, the stationary probability 
distribution density is often unknown and cannot be used as preknowledge during the simulation. 
We compare the algorithms efficiency with standard Brownian dynamics simulations and other WE 
methods. 
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Rare events are ubiquitous in many biological, chemi- 
cal and physical processes [l|, ■ Whereas the density of 
states is known in systems at thermal equilibrium, inter- 
esting phenomena often occur in non-equilibrium systems 
Unfortunately, many such problems are inaccessi- 
ble to analytic methods. Therefore computer simulations 
are a widely used tool to estimate the density of states 
or transition rates between them 4 , [j| . Since standard 
Brownian dynamic simulation 0, 7| provides computa- 
tional costs that are inversely proportional to the state's 
probability, specialized methods [8l-[lQ| have to be used 
to adequately sample rare events, i.e. states with low 
probability or low transition rates. 

In the last decades, flat histogram algorithms have 
been developed, allowing one to evenly sample states with 
highly different probabilities. These algorithms are im- 
plementations of the umbrella sampling [l2| , where each 
state is sampled according to a given probability distri- 
bution, the so-called umbrella distribution. Within non- 
equilibrium umbrella sampling (NEUS) [l3| the space of 
interest is divided into different but almost evenly sam- 
pled subregions. The interaction between different re- 
gions occurs solely due to probability currents between 
then. Whereby the probability distribution within a re- 
gion is then calculated by performing Monte Carlo sim- 
ulations. 

In order to calculate low rates between a starting and 
a final state, forward flux sampling methods can be used 
(for a review see) [lij]. These methods introduce a se- 
quence of surfaces between these states and introduce 
walkers (copies of the system) to perform weighted trajec- 
tories according to the underlying dynamics. If walkers 
cross one of the surfaces, getting closer to the final state, 
new walkers with smaller weights are introduced. Finally, 
many walkers with particular small weights reach the fi- 
nal state. The consideration of the particular weights 
allows one to calculate very low rates in a finite simu- 



lation time. Recently, extensions to these methods have 
been developed to calculate both, transition rates using 
umbrella sampling (l5j and probability distributions us- 
ing forward flux sampling |16j algorithms. 

In this work, we present an algorithm, based on the 
previously developed weighted-ensemble (WE) Brownian 
dynamics simulations jl7H20| . that allows one to calcu- 
late the stationary probability density function (SPDF) 
as well as transition rates between particular states. Like 
in WE simulations the space of interest is divided into 
several subregions and the probability for finding the sys- 
tem in them is calculated by generating equally weighted 
walkers in each region. By moving to the underlying dy- 
namics, the walkers transport probability between the 
subregions. Thus, WE methods are usually applied to 
systems of Brownian particles moving in a potential land- 
scape dum. 

We are interested in an algorithm which allows sim- 
ulations of stochastic dynamical systems which, apri- 
ori do not obey detailed balance for probability fluxes 
or suppose some special topology of the flow [22j - [26j . 
Such are given for Brownian particles in conservative 
force fields under the influence of additive noise. Even 
canonic dissipative systems possess a vanishing probabil- 
ity flow if transformed to the energy as dynamic variable 
[24 l27l [28j . In consequence, both systems allow exact 
analytic solutions of the SPDF. In contrary, we aim to 
develop an algorithm which does not assume that neither 
the deterministic nor the stochastic items (see Eq. (Q]) be- 
low) undcrly such conditions. Thus, no information on 
the SPDF can be used for the simulations. 

In general the algorithm can be applied to arbitrary 
dynamical systems of the form: 



x n = /„(x) + 5„(x)£„(t), n = 1, d, 



(1) 
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where d is the number of stochastic time-dependent 
degrees of freedom x n {t),n = l,...,d; the func- 
tions /n(x) describe the deterministic velocities 
for the n-th direction; £n(i) represents zero-mean 
Gaussian noise with delta-like correlation function 
(£.n(t)£,m(t')) = 5 nm 5(t — t'). The noise intensity along 
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Figure 1. Scheme of the interpolation points xi and the sup- 
porting points Xj used in the numerical calculations. In this 
scheme we have used [ A j 1 ] = 4. 



the n-th direction is scaled by the functions g n (x) , which 
in general depend on the vector x — (xi, X2, Xd) T - 
We are interested in high precision sampling of the 
stationary probability current J s t(x) and the SPDF 
p s t(x) of finding the system in the d-dimensional cube 
[x\,xi + dx\\, [xd,Xd + dxd] with a finite resolution. 
We will specify the resolution by the number M n ^ res of 
evenly spaced supporting points along the n-th direction, 
for which we will determine p st (x) . 

This article is organized as follows: In section U we in- 
troduce an algorithm, based on WE methods [17] . that 
allows one to calculate low probabilities and rates. After- 
wards, we study one- and two dimensional model systems 
and analyze the algorithms efficiency compared to Brow- 
nian dynamics simulation (BDS) and WE techniques. 



I. THE ALGORITHM 

First, for sake of clarity, we restrict ourselves here 
to particle motion in one dimension, d = 1, under ad- 
ditive noise. The deterministic part of the dynamics 
Eq. ([T]) can be always represented as a conservative force 
f(x) — —U'(x). The noise strength is scaled by the pa- 
rameter D, we put g(x) = \2D. For such systems, the 
SPDF is known to be 

Pat (x) = Z~ t x e"^, (2) 

where Z st is a normalization constant. 

We are interested in finding numerically the system's 
SPDF, p s t(Xj), at a set of M res evenly spaced supporting 
points Xj in a finite part of the physical space, given 
by x G [L~,L + [. The region of interest is divided into 
M > M res subregions of size Ax = L+ jU L , the i-th 
subregion is bounded by (xj, Xi+i), i = 0, . . . , M— 1, with 
Xi = iAx + L~ . Supporting points are given explicitly 
by Xj = L~ + (j - |)AI res , j = 1,2,..., M res , with 
AX res = ~\Ax, see Fig. Q] Here \z\ denotes the 

largest integer smaller than or equal to z. 

Let us introduce the probability Pi(t) for finding a par- 
ticle in the i-th subregion at time t, and the correspond- 
ing set P(i) = (P (i),Pi(f), P M -i{t)). Initially, no in- 
formation on the system is available, thus, each subregion 
is given an arbitrary amount of probability Pi(0), simply 
fulfilling the normalization condition J2i=o 1 = 1- 

Naturally, equilibration can be accelerated if one already 
has information on the SPDF of the system (Eq. j!}). 



In that case, one can choose P(0) close to the set P st , 
which optimally approximates the SPDF: 

r x i+i 

Pst,i = dx p st (x). (3) 

J Xi 

However, in general no such information is required. 

A. Time evolution 

After setting the initial set P(t = 0), the time evolution 
of the P(t) — ^ P(< + h) is performed using three different 
steps. 

We start with a redistribution step, in which N walkers 
(copies of the system) are uniformly distributed in each 
subregion. Besides their individual positions x^(t), where 
i = 0,...,M — 1 denotes the particular subregion and 
k = 1, N the individual walkers, each walker possesses 
a given amount of weight q%{t). This is nothing but the 
present probability in the i-th subregion distributed to 
the N walkers, which yields 

flf (*) - (4) 

Note that one does not need to introduce walkers in sub- 
regions with Pi(t) = 0. 

After the redistribution step has been performed, Eq. 
([T]) is integrated for all walkers, using a Brownian dy- 
namic simulation step h and an arbitrary integration 
scheme. This integration step realizes the time evolution 
x\(t) — > x\(t + h). Here walkers transport probability 
between the subregions. As walkers are independent of 
each other, it is of importance to note that the particular 
time evolution of each one of the N x M walkers is due 
to different sample paths in the stochastic parts of the 
Langevin equation. 

Lastly, an updating step is performed, in which the new 
probabilities Pi(t) — > Pi(t+h) are calculated by summing 
up the weights of all walkers that are currently located 
in the particular subregion, 

m+h)= J2 «*(*)• ( 5 ) 

i' ,k\x k l {t+h)^ 1 {xi,Xi^\) 

In what follows, we will name the sequence of redis- 
tribution, integration and updating step as running step. 
After an equilibration time Tthenm the set P(t) reaches 
a stationary regime, where the Pi(t)'s fluctuate around 
their mean values (Pi). 

B. Calculating the stationary probability 

The individual (Pj) are estimated by averaging over 
a total amount of Nt sets P(t^), t = 1, 2, Nt, taken, 
after the system has reached the stationary regime, every 
n av running steps: tt = T t herm + {£ — l)n av h. It turns 
out that the mean probabilities (Pi) coincide with the 
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stationary probability P st .i (see Eq. (J3J), for compatibly 
chosen time step h and the size of a subregion Aa; (see 
Sec. EH). 

Finally the SPDF on the supporting points p s t(Xj) is 
calculated by adding the adjacent (Pi) and dividing by 
the size AX res (in order to have a properly normalized 
PDF): 



Pst{Xj) 



1 



AX r 



E 

=o-i)ri^ri 



(6) 



C. Calculation of the probability current 

The stationary probability current J s t(x) at position 
x can be easily calculated by adding up (with the right 
sign) the weights of all walkers, passing x to the right 
and to the left per unit time. In practice x should be the 
boundary of a subregion. If x = x% the current J(x,)(i) 
is given by: 



j(xm = U E - E «?'(*) 



(7) 



i'k&Li 



and i?i indicate these walkers which cross the boundary 
moving rightwards, i.e. x^it) > Xi A xf,(t — h) < x%. 
Alternatively, Li assign walkers transporting weight left- 
wards, x^, (t) < Xi A x\, (t - h) > Xi. 

Averaging over Nt such estimates, taken in the sta- 
tionary regime, leads to the average current (J(xj)), 
which converges towards J 8 t(xi) for Nt — > oo. 



D. Implementation of boundary conditions 

The implementation of boundary conditions for the 
probability current or the SPDF is straightforward. 
Right now, absorbing boundaries are already imple- 
mented at L~ and L + , since walkers that pass these 
boundaries are not located in any subregion. Therefore, 
their weights will get lost in the next updating step. Re- 
flecting boundary conditions at L + can be implemented 
by setting x*{t + h) -> 2L + — x$(t + h) for all walkers 
with x\{t + h) > L + . Hence, the probability current 
at L + will be zero. Reflecting boundaries at L~ can be 
implemented analogously. 



E. Convergence criteria 

1. One- dimensional systems 

In order to ensure, that (Pi) and (J(xj)) converge 
towards the stationary probability distribution and the 
probability current of Eq. (JTJ) for Nt — > oo, the time 
step h and the size of a subregion Ax have to fulfill 



specific criteria. This is due to the redistribution step, 
where walkers are uniformly distributed in each subre- 
gion. This implicates statistical errors, since they can 
reach positions in a subregion that, are inaccessible or 
at least more improbable. Therefore, walkers can more 
easily escape from potential minimums or reach regions 
of low probability. This effectively flats the probabil- 
ity distribution, leading to more probability in regions of 
low probability, for instance, around local maximums of 
U(x), and less probability in the potentials minimums. 
In order to overcome this problem, earlier works [13, EH 
have stored the positions and weights of all walkers. In 
the next redistribution step, walkers were only spaced on 
the stored positions according to the weights belonging to 
them. This requires a lot of computer memory, especially 
for large N and M. 

However, we found that one does not need to store 
these information, if the subregions are small enough to 
ensure that walkers have a non-negligible probability to 
leave them during one integration step. As a measure 
of how far a walker can step, due to the fluctuations, in 



one time step, we use the diffusion length Ldif = 2 



/ 

Thus, the size of a subregion Ax should be small com- 
pared to the diffusion length L^f 



Ax < 2\ 



(8) 



The distance a walker can pass during an integration step 
is not only determined by L^if, but also by the determin- 
istic dynamics, leading to a step length L^et = f{x)h in 
first order. Usually /(x) changes very fast at the bound- 
aries of the simulation area, which produces high deter- 
ministic velocities and regions of low probability. Walkers 
can only reach these regions, if the fluctuation arc strong 
enough to balance the deterministic force, i.e. if 



\f(x)\h<2y 



(9) 



for all x G [L ,L + [. This leads to a condition for the 
time step h: 



h < hr, 



AD 



max l£ [L-,L+] / 2 (x) 



(10) 



Hence, lower time steps allow one to sample regions, far 
from the potential extrema and for instance, the tails of 
the SPDF. If a larger time step is chosen, the (Pi) will 
run to zero in subregions with larger deterministic force. 

Since it is often difficult to fulfill Eq. (fTOj) in the entire 
simulation area, one should choose a time step, which 
allows one to fulfill Eq. ©. 



2. Multidimensional systems 

In general, our method can be applied to stochastic dy- 
namical systems in arbitrary dimension d > 1 (Eq. (|T|)). 
For the foundation of the algorithm, we refer to the Ap- 
pendix Sec. [Bj However, in order to ensure convergence 
in a finite region A, the criteria for the time step h and 
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the size of a subregion Axi along the i-th direction should 
be fulfilled. For noise dominated directions the criteria 
hold, therefore, the condition for the time step (Eq. ([TO]) ) 
becomes: 



h m min ( 



4£>„ 



i<n<d max j4 (/„(x)) 5 



(11) 



and the criteria for the size of a subregion along the n-th 
directions (Eq. (|HJ|) reads: 



Ax n < 2^/D n h. 



(12) 



However, there might be directions without any noise 
(D n =0). In that case the length of a subregion should 
ensure, that walkers can leave it due to the deterministic 
term /„(x), leading to 



Ax n < f n (x)h. 



(13) 



otherwise information on the deterministic dynamics gets 
lost during the redistribution step, since walker that stay 
into a subregion do not produce any change in the Pi 
and are again randomly placed in their subregion during 
the next redistribution step. For equally sized subregions 
Ax n should be the minimum value of Eq. (|13|) with re- 
spect to all x in the simulation area for each of the d 
directions. These criteria can lead to a huge number of 
subregions, especially in high dimensional spaces. 

To appropriately reduce the number of subregions, the 
Ax n should be chosen in order to locally fulfill the crite- 
ria. This was implemented by a grouping algorithm, that 
groups original subregions into "larger" ones as long as 
walkers can leave these due to the deterministic term (Eq. 
1131) . Depending on the system, this procedure highly re- 
duces the total number of "larger" subregions M group . If 
the grouping algorithm was used, N walkers are randomly 
placed in each of these "larger" subregions and a proba- 
bility Pi of finding a walker in the corresponding area was 
introduced. We find that such grouping highly reduces 
the computational costs, since less subregions and there- 
fore less walkers are required. Since walkers jump out of 
these regions until the next redistribution step starts, the 
algorithm still approximates the correct SPDF. 



F. Simulation techniques 

Simulations were performed on a Intel®Xeon ®CPU 
E31245 @ 3.30GHz processor with 16 Gb DDR-3 RAM. 
The algorithm described above was implemented in a 
C ++ program for one- and two-dimensional systems. 
Runs of the algorithm are specified by the time step h, 
the size of a subregion Ax, (Ay, in two-dimensional prob- 
lems), the simulation area, given by L~ and L + (L^, 
Ly), the number of walkers per subregion N, the ther- 
malization time Ttherm-, and the number of running steps 
between two sets of P denoted The numerical 

integration of the Langevin Eq. (JTJ) was done using a 
Heun scheme. A resolution of M res = 200 was used in 
any direction. 



To compare the results with other methods, we also 
perform Brownian dynamics simulation (BDS) using 
Nsrown initially uniformly placed particles in the simu- 
lation area. Integration was done using the Heun scheme 
0, HH with integration time step h Brown- After a ther- 
malization time T t herm, Brown the particles positions were 
recorded after time intervals At Brown- The BDS was 
given a running time T run (real CPU time) which usu- 
ally equals the time our algorithm needs to produce its 
results. After T run the BDS was stopped and the SPDF 
was calculated using the recorded particle positions. We 



setT, 



therm, Brown 



Tti 



therm-i h Brown 



— /i, AtBrown ^qii^ 



to make results comparable. 



II. MODEL SYSTEMS AND RESULTS 

A. One dimensional system 

In order to demonstrate the implementation of the al- 
gorithm, wc study ovcrdamped Brownian motion in a 
bistable potential U(x) = — ^- + Correspondingly, 
we put f(x) — x — x 3 and g(x) = \J1D in Eq. (jTJ) which 
results in a bistable system which is often used to s tudy 
bistable systems or stochastic resonance therein pll. l3Cll]. 
The two stable states come up to the potentials mini- 
mums, located at x = — 1 and x = 1, respectively. The 
corresponding SPDF is given by Eq. ([2]). For low noise 
strength, the SPDF attains sharp peaks at the potentials 
minimums and decreases down to low values at the bor- 
ders and the local maximum, for instance for D = 0.01 
Pst (0) « 10- 11 . 



1. Equilibration 

At first, we study the equilibration process, performing 
simulations with different numbers of walkers per subre- 
gion N . Results are shown in Fig. [2] Analyzing the 
time dependence of the probability Pi(t), we find that 
longest thermalization time occurs at the local maximum 
of U{x) . Note that runs with larger N thermalize at lower 
t, but one needs more integration steps. 

After thermalization has been achieved, we evaluate 
the coefficient of variation of the probability, given by 



c(xi) 



(Pi) 



(14) 



where averages (■ • • ) are performed for a fixed number 
Nt of sots P and different N and Xi. Wc find it to scale 
according to (data not shown). 



2. Stationary probability density function 

We start to calculate the SPDF p st in the region 
[L~,L + [ for a small noise strength (D = 0.01). The 
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Figure 2. PDF at x = obtained from Pi by PDF(t) = 
for the subregion containing x — for D = 0.01 and different 
numbers of walkers per subregion N. Parameters are chosen 
as in run 1 (see Tab. U). 
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Figure 3. Estimates for the stationary probability density 
obtained from the algorithm for run 3 (see Tab. |TJ) and by 
using a Brownian dynamics simulation for D — 0.01. Analytic 
results are obtained from Eq. l[2]). 



run 


L~ 


L+ 


h 


Ax 


M 


1 


-1.4 


1.4 


0.011 


0.00104869 


2670 


2 


-1.75 


1.75 


0.0015 


0.000387297 


9037 


3 


-2.5 


2.5 


0.0001 


0.00010775 


46404 



Table I. Time step h and box size Ax according to convergence 
criteria Eq. (|10|) and where we choose h = fe, y x and 



Ax = iL d 



time step h and the box size Ax are set according to the 
criteria (see Tab. fl}. Using the results shown in Fig. [5J 
we set the thermalization time T t herm = 50 for a run 
with N = 2. Time averages after thermalization were 
calculated over an ensemble of Nt = 10 4 sets P. Results 
for the SPDF are shown in Fig. [3J The algorithm cal- 
culates the tails of the distribution down to 10~ 300 cor- 
rectly, after a running time T run « 27 hours. We also 
calculate the SPDF using Brownian dynamics simulation 
using N Brown — 10 4 , which stops estimating at a level of 
10~ 6 after the same running time. Further runs were 
performed (see Tab. |T](run 2) and (run 3)), approximat- 
ing the tails down to 10" 10 (M group = 1136) and 10~ 48 
(M group = 6789) after a running time of w 30 sec and 
sa 10 min , respectively (data not shown). 



Simulations for different values of 



indicate that 



insignificant deviations from the analytic SPDF occur for 
Ax > 4t- 



tions [3l| and in the field of neuroscience [32J . 

Initially, only N particles are assigned at the subregion 
including x m i n — 1, so approximating an initial delta- 
like probability distribution for t = 0. Furthermore, 
an absorbing boundary right behind the local maximum 
(x a bs = —0.01) is included. To fulfill normalization of 
the SPDF, walkers that reach x a bs are reinjected imme- 
diately at x min . The escape rate to pass the barrier is 
given by the probability current J(x max )- For small noise 
intensities, the probability current on top of the potential 
barrier can be described using Arrhcnius law, namely: 



J(x ) oc e 



(15) 



where AU = U(x max ) — U(x m i n ) — 0.25. Since strong 
fluctuations are rare, but possible, we will use J{x a b s ) to 
approximate J(x max ). Probability currents J(x a b s ) were 
recorded for each time step and averaged over a sequence 
of n av = [%-] running steps, resulting in (J(x a b s )}- 
Figure|4]shows the time dependence of ln( J(x a bs)}- After 
a relaxation regime, where the current decays exponen- 
tially, the current reaches its stationary value. The values 
of ln( J(x a bs))i averaged over the stationary regime, are 
shown in Fig. [5] for different noise intensities. Fulfilling 
the criteria described above, the algorithm reproduces 
well Arrhenius law down to hi(J(x max )) ~ —650 corre- 
sponding to a current J{x max ) ~ 10 -286 . 



Efficiency compared to weighted-ensemble Brownian 
dynamics simulation 



3. Probability current 

Next, we present that our algorithm can be used to 
calculate the escape rate to pass the energy barrier at 
Xma X = 0. Such problems are typical for chemical reac- 



In order to compare the efficiency of two algorithms, 
important quantities are the transient time required to 
first reach the steady state. Once the algorithm reaches 
the steady state, we quantify the size of the fluctuations 
by the coefficient of variation Eq. (jT4")) at the potential's 



6 



o 

-100 
-200 
-300 
-400 
-500 
-600 
-700 



HI 1 " *!!' ** ! 11 



D=0.1 

D=0.004 

D 0.002 

D 0.0008 

D=0.0005 

D=0.0004 ; 

100 200 



10 



300 
t 



400 



500 



600 



Figure 4. Time dependence of the probability current J(x a t B ) 
obtained from our algorithm for L~ — —0.02, L + — 1.39 and 
decreasing noise intensities (from top to bottom). The time 



Ldif- Note that h„ 



step is set h = "™ a " and Aa 

and L,uf vary according to Eq. (|iop and respectively, 
resulting in larger running times for smaller D. 
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Figure 5. Stationary probability current as a function of the 
inverse noise strength obtained by time averaging the data 
partly shown in Fig. [4] in the stationary regime. Error bars 
show three standard deviations of the stationary data. 



local maximum x = 0. The computation time mainly 
depends on the number of integrations Ni nt needed, to 
reach the stationary regime. In order to compare the ef- 
ficiency, the size of fluctuations (Eq. (fH)) ) in the local 
minimum relative to p st (0) w 3.88717 x 1CT 11 during the 
stationary regime is plotted over Ni nt in Fig. [S] The 
most effective algorithm would be located close to the 
origin. Comparing the efficiency of standard WE simu- 
lations and our algorithm, we find that WE simulations 
with low N equilibrate faster. However, the precision 
highly depends on the fluctuations during the stationary 
regime. To produce results of same precision (same c(0)) 
both algorithms approximately need the same Ni nt . Usu- 
ally WE simulations were performed using thousands of 
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Figure 6. Relative fluctuations plotted over the number of 
integration steps, needed for equilibration, for the algorithm 
(black) and standard (WE) simulation (red). The number of 
walkers per subregion is shown for each point. Simulations 
were done for run 1 and 2 (compare Tab. [T|. 



walkers per subrcgions, resulting in a low value of c(0). 
Here runs of our algorithm producing the same c(0) need 
much less N and have memory requirements independent 
ofN. 



B. Two-dimensional systems 

1. Poincare Oscillator 

As an example of a two-dimensional system with 
known SPDF, we consider the Poincare oscillator [27], 
HI, HI| , represented by the dynamical system: 
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v 

(a 



2\ 

v )y 



(16) 



where represent delta correlated white Gaussian 
noise with zero mean. Using the energy function 
H(x, y) — \(x 2 +y 2 ), which only depends on the distance 
to the origin, one can calculate the associated SPDF: 



Pst(x,y) 



cxp( 



r aH (x,y) - H 2 (x,y) , 
D 



(17) 



(see appendix Sec.[X|. Since noise only applies to the y- 
direction, the lengths of a subregion Aa; and Ay in x- 
and y-direction are calculated by Eq. (fT3|) and ©, re- 
spectively. The minimum of \ f x {x,y)\ = \y\ is equal to 0, 
therefore, we choose Ax = Ay h, which corresponds to 
a first order approximation of \f x (x,y)\ in the next sub- 
region. Results for the SPDF are shown in Fig. [71 In- 
terestingly, the algorithm approximates better the SPDF 
along the direction where no noise was applied. Here the 
SPDF is sampled down to 10~ 30 . We found that the al- 
gorithm slightly oversamples the analytic SPDF in the 
tails. This is due to the statistical errors, made dur- 
ing the redistribution step. By reducing the size of a 
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Figure 7. Contour plots of the SPDF (D = 0.1 and a = 1) 
obtained from the algorithm (black), using the parameters 
L x = L y — —3, = Ly = 3, Ttherm = 8, iVr = 10 3 , h = 
0.01, M a = 189737, M y = 1897, M gro »p = 1688940, AT = 2, 
(Ay = ^Ldif) and the analytic solution (red) Eq. (|17[) . Con- 
tour lines are labeled according to represented values of the 
SPDF and show the rotational symmetry. The SPDF pos- 
sesses its global maximum at H(x,y) — ^, corresponding to 
the unit circle for our choice of a, and a local minimum in the 
origin. Running time ~ 5 hours. 



subrcgion, this error can be reduced further. Along the 
y-direction, noise is applied. Here the behavior is similar 
as in the one dimensional example (see above). For runs 
with larger Ay (results not shown) the algorithm slightly 
oversamples the SPDF in the origin. 



2. Bistable system with colored noise 

As a further example, we calculate the SPDF of the 
two dimensional system: 



x = x — x + y 



where r denotes the time scale separation between x and 
the colored noise y. The white Gaussian noise has 
been already described above. This system has been 
studied previously in [HI, |3(| • Like in the bistable sys- 
tem we have studied above, the SPDF has maximums at 
(x, y) = (1,0) and (x,y) = (—1,0). However, for some 
combinations of r and D, the SPDF possesses a local 
minimum at (x,y) = (0,0). Plots of the SPDF are de- 
picted in Fig. [SJ Our algorithm samples the SPDF down 
to 10 -12 , whereas BDS breaks down at a level of 10 -6 . 



Si 




-0.5 



-1.5 L 



-1.0 -0.5 0.0 0.5 1.0 

X 



Figure 8. Contour plots of the SPDF of Eq. fig]) (D = 0.1, 
t — 2 and r = 2.5) obtained from the algorithm (black) 
and BDS (red) of the same running time. Parameters: (top) 



— Ly 

0.0889, M x 



14319, M v = 1273, M„ 



-1.5, Lt = L+ = 1.5, Ttherm = 10, N T = 10 3 , 

= 165003, N = 
10 3 , and (bottom) 
1.5, T therm = 12, N T = 10 3 , 
h = 0.0889, M x = 17899, M y = 1591, M group = 206957, 
(18) TV = 2,(Ax = Ay h, Ay = ±L dif ), 7V Sro „„ = 10 3 . Running 



2,(Az = Ay h, Ay = ^L di; ), N Br 
L~ = L„, = 1.5, Lt~ = L 



times are 19 min (top) and 21 min (bottom) 



The minimum, occurring for r = 2.5 was clearly found 
by the algorithm. 



3. FitzHugh-Nagumo-system 

As a last example, we consider the widely used 
FitzHugh-Nagumo-system [37} , which is often used in 



the field of Neuroscience [38[ or to study synchroniza- 
tion [3!l[40[ and coherence phenomena 41 1, represented 
by: 



D = 0.01 



x 3 — y) 



y 



jx-y + b+^2L\^(t). 



(19) 



Here e denotes the timescale separation between the acti- 
vator variable x and the inhibitor variable y. £x(t), £y(t), 
represent independent zero-mean delta-correlated Gaus- 
sian white noises. We want to study the stationary prob- 
ability density in the case of D x = D y = D for a time 
scale separation e = 0.1. We set the parameters ac- 
cording to Ref. H3 to b = 1.4, e = 0.1 and 7 = 2. 
Thus, the system is in the excitable regime. Since the 
deterministic part of the equation for the activator vari- 
able increases very fast if x is increased, we have to 
choose a time step h = 0.01, which is small enough, that 
the walkers' steps are small compared to 1, but allows 
us to fulfill the criteria for the size of the subregions. 
Contour plots of the SPDF are shown in Fig. [9] Es- 
pecially regions of low probability are much better sam- 
pled, using the algorithm. In the case of low diffusion 
(D = 0.01, Fig. [9] (top)) the algorithm runs down to 
10~ 16 , whereas BDS stops at a level of 10 -6 . Especially 
the minimum is much better sampled by our algorithm. 
Note that the local maximum located in the surroundings 
of (x,y) = (0.7,0.5) was not found by BDS. For higher 
diffusion values (D = 0.1, Fig. |9] (top)), significant dif- 
ferences can be found only in the tails. 



III. DISCUSSION 



A. One dimensional system 



In the one dimensional system (see section III A[) we 
succeeded to approximate the SPDF down to 10 -300 . 
If equally sized subregions are used, a huge number of 
subregions has to be implemented, in order to fulfill the 
convergence criteria (see section ITEf . However, by eval- 
uating the SPDF according to Eq. (j6|) the additional 
computational costs also reduce the fluctuations of the 
estimated SPDF. If the subregions are to large, the the- 
oretical SPDF is overestimated by the algorithm in po- 
tential minimums and underestimated in the potentials 
maximums. 

We also managed to calculate escape rates of size 
10 -286 . Although the finite size of a subregion leads to 
small errors, in the calculated SPDF, Arrhenius law is 
well reproduced for such small probability currents. 

At the stationary regime the Pi fluctuate around their 
mean value, estimating the SPDF. By either increasing 
the number of averages Nt or the number of walkers per 
subregion N these fluctuations can be reduced, leading to 
higher precision. The estimation error scales with V 'Nt 
and \/N, respectively. However, increasing the number 
of walkers per subregion also effects the computational 
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X 

Figure 9. Contour plots of the the SPDF of Eq. fod 
(D = 0.01 (top) and D = 0.1 (bottom)) obtained from 
the algorithm (black) and BDS (red) with the same running 



time. Parameters: (top) L x = L y 



"2, L x — Li — 2, 



N 



= 5, N T = 10 d , h = 0.01, M x = 4000, M y = 4000, 



— 9 f Ax = I) 



(bottom) L x 



Nsrown — 10' 
= ~ 2, L x = Ly 

1333, M v 



Mgroup = 314027, and 



2, T, therm ^? 

1333, N = 2, 



M group = 146568. 



N T = 10 J , h = 0.01, M : 

(Ax = Ay = jfiLdif), N B rown = 10 

Running times are 27 min (top) and 19 min (bottom). Sim 
ulation times for the BDS were chosen three times larger. 



costs during the thermalization. Simulations for differ- 
ent N show that runs with higher N become stationary 
faster, but this does not compensate for the additional 
computational costs. We also find, that increasing N 
slightly improves the sampling of the SPDF's tails. Once 
the system is in the stationary regime the increase of 
the computational costs scale linearly with N and Nt- 
An advantage of simulations with small N is, that one 
does not need to perform running steps for subregions 
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with Pi — and the number of such regions naturally in- 
creases for small N. We usually use the smallest possible 
N = 2 and scale the estimation error by increasing Nt- 

B. Comparison with weighted-ensemble Brownian 
dynamics 

Since the general idea of our method was adapted from 
prior simulation tec hniq ues known as weighted-ensemble 
(WE) sampling [13, Il8l . l20j , we want to discuss advan- 
tages and disadvantages of our algorithm in comparison 
with these techniques. The main difference in WE tech- 
niques is the redistribution step. Here, using WE tech- 
niques, positions and weights of all walkers are stored, 
and new walkers are introduced on the stored positions 
considering their particular weights. 

In our algorithm, there is no need to store any position 
or weight, since walkers are randomly placed in each sub- 
region. The resulting statistical errors can be neglected, 
if the size of the subregions fulfills conditions which, un- 
fortunately, lead to much larger numbers of subregions. 
By averaging the probabilities, these extra computational 
costs contribute to an reduction of fluctuations in the sta- 
tionary regime. 

The comparison of the computational costs until equili- 
bration and the achieved precision shows, that both algo- 
rithms posses the same efficiency for high precision runs. 



the diffusion length. In contrast to WE methods, the re- 
quired memory does not depend on the number of walk- 
ers, which leads to less memory requirements for runs 
with large numbers of walkers. 

Special attention was payed to non-equilibrium dy- 
namical systems. Applying the method to one- and two- 
dimensional model systems, we analyze its efficiency com- 
pared to standard Brownian dynamics simulation. Our 
method outperforms Brownian dynamics simulation by 
several orders of magnitude and its efficiency is compara- 
ble to weighted-ensemble Brownian dynamic simulations 
in all studied systems and lead to impressive results in 
regions of low probability and small rates. 
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C. Two dimensional systems 

In the case of two dimensional systems, we found that 
the algorithm outperforms Brownian dynamics simula- 
tions. However, the number of subregions needed ac- 
cording to the criteria can be really high, especially if 
there are some directions without any noise. Here it is 
necessary to size the subregions to fulfill the criteria lo- 
cally. A first step in that direction has already been done 
by the implementation of the grouping algorithm. We 
are quite confident that it is possible to reduce the com- 
putational costs by optimizing the Mesh. The analyzed 
examples show, that the running times highly depend on 
the investigated system. For the bistable system with 
colored noise and the FitzHugh-Nagumo-system, we ob- 
tained good results within « 20 minutes, whereas the 
algorithm needed about 5 hours for the Van-der-Pol os- 
cillator. 



IV. CONCLUSION 

We provided and tested an algorithm that allows the 
calculation of low probabilities and low rates. The al- 
gorithm is based on WE Brownian dynamic simulations, 
but uses a uniform distribution of walkers within each 
subregion. To our findings, the resulting statistical errors 
can be neglected if one uses subregions small compared to 



Appendix A: Stationary probability density of the 
Poincare oscillator 



Using the energy function H(x, y) = \{x 2 +y 2 ), we get 
from Eq. (|I6[) to a representation as a canonic dissipativc 
system: 



x 

y 



dyH 



d y (aH - H 2 ) - d x H ■ 



(Al) 



The corresponding Fokker Planck equation in the x, y 
phase space for the SPDF p s t(x, y) reads (271 ]: 

dtPst = = -d y Hd x p st + d x Hd y p st 

-d y (d y (aH -H 2 ) Pst ) + Dd 2 Pst . { ' 

Using the ansatz p st (x,y) = Ps t(H(x,y)), the first two 
items at the r.h.s. cancel. The remaining second line can 
be integrated once. Assuming an exponentially decaying 
SPDF at infinitely large energies yield the disappearance 
of the irreversible probability flux in y-direction (22), HH . 
One finds, afterwards : 

P«±{aH-ll*)=D±p«. (A3) 
It leads to Eq. (fTTl) . 
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Appendix B: Foundation of the algorithm 

In this Appendix we show that the presented algorithm 
is described by a corresponding Master equation for the 
probability distribution density Pi{t) for the case of a 
Markovian hopping process between boxes. We use a sin- 
gle index i to label the boxes, but the argument applies 
to any spatial dimension d. We identify the dynamics 
of the stochastic system which shall be simulated with 
the discrete stochastic dynamics of the walkers. The lat- 
ter is defined via the matrices of probabilities per unit 
time w(i — > i') which describe the hopping in the given 
discretized space. We assume that it shall converge for 
sufficiently small time scales and box lengths to the out- 
going dynamics. 

Let assume that we have the probability distribution 
density given at time t in every box with index i. It holds 

= 1. (Bl) 

i 

We redistribute the probability in every box to N walk- 
ers. In result any walker k = 1, . . . , N of the same box 
gets an identical weight 

rf(*) = ^- (B2) 

until the next new redistribution. 

At later time t + h the walkers may be still inside the 
box or may have jumped to other boxes. Let Ui denote 
all possible box-indices which can be reached during a 
single step from the i-box. Then in accordance with the 
assumption above, the probability is w(i — > i')h that 
during h a single walker leaves i-box and jumps to the 
box with index i' € Ui. By this hopping the walker 
k carries the weight q^(t) to the new box, which step 
is, obviously, connected with a lost in the outgoing box. 
Therefore, the lost per particle for this specific hopping 
from i — > i' can be expressed as 

w{i -^i')hq^(t). (B3) 



The whole lost of weight will be realized on all possible 
hopping channels. It is identical for all N particles being 
located in the present box. Hence, the full lost becomes 

Nq*{t) w(i -> i')h = Pi(t) ^ w ( { -*i')h. (B4) 

i'eUi i'eUi 

Alternatively, one can also introduce U[ as the boxes 
where from walkers can reach the i-box. Then, gain of 
weight transferred by every walker arriving at the i-box 
is expressed, if i' e U-, by 

w{i'^i)hq^, (B5) 

Again, summing over the different hopping steps and con- 
sidering that all the TV walkers inside the box with i' E U[ 
reach the i box, yields the gain 

-> i)hNq$(t) = ^ w(i' -> i)hP v {t). (B6) 

Therefore, the balance of transported weight results in 
the following shift of the full weight in the i-box at time 
t + /;, compared to the former one 

Pi(t + h) - Pi(t) = - Pi(t) w ( l "> *') h 

+ w(i' ^i)hP t ,{t), (B7) 

plus corrections of order 0(h 2 ) corresponding to cases in 
which two or more walkers jump from one box to an- 
other during the time interval h. After dividing by the 
time step h and taking the limit h — > we obtain the 
wanted Master equation for the discretized dynamics in 
the boxed phase space 

d t Pi(t) = - -> *") 

+ J2 «>(*'-> *)iV(*)- (B8) 
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